Stochastic Resetting for Enhanced Sampling

We present a method for enhanced sampling of molecular dynamics simulations using stochastic resetting. Various phenomena, ranging from crystal nucleation to protein folding, occur on time scales that are unreachable in standard simulations. They are often characterized by broad transition time distributions, in which extremely slow events have a non-negligible probability. Stochastic resetting, i.e., restarting simulations at random times, was recently shown to significantly expedite processes that follow such distributions. Here, we employ resetting for enhanced sampling of molecular simulations for the first time. We show that it accelerates long time scale processes by up to an order of magnitude in examples ranging from simple models to a molecular system. Most importantly, we recover the mean transition time without resetting, which is typically too long to be sampled directly, from accelerated simulations at a single restart rate. Stochastic resetting can be used as a standalone method or combined with other sampling algorithms to further accelerate simulations.

parameter of 100 ∆t. The condition for the FPT was checked every 100 time steps. We used the AMBER99SB force field. Available GROMACS input files for this system S3 were converted to LAMMPS format using Intermol. S4 With them, we obtain a free energy difference between the conformers of ∼ 9 kJ mol −1 , which is in reasonable agreement with reference values. S5 The FPT distributions with and without resetting were obtained from 10000 trajectories each. For the trajectories without resetting, 307 trajectories did not show a transition within 4 µs. They are not plotted in the probability density of Fig. 1 but are included in calculating the mean FPT and speedup. We note that, as a result, the speedup gained by resetting that we report for alanine dipeptide is a lower bound.

Implementation of stochastic resetting
Stochastic resetting (SR) was implemented in the input files as explained below. A stopping mechanism after the first transition was also incorporated through the LAMMPS input. The initial velocities, and their values after each reset event, were sampled from the Maxwell-Boltzmann distribution at the relevant temperatures using Python. For Poisson resetting, waiting times between resets were also sampled using Python, from an exponential distribution, f (τ ) = re −rτ , where r is the restart rate. Below, we give a simplified example of the implementation for a two-dimensional simulation with only three reset events for clarity.
The initial position in this example is fixed at (1, 0)Å and the first transition (passage)

Laplace transforms for the inverse Gaussian distribution
We used Eq. 1 for the inverse Gaussian probability density function. The expression for its analytical Laplace transform is given in Eq. 2. We used the parameters L = 1000, V = 1 ps −1 and D = 12500 ps −1 , which lead to a mean FPT of 1 ns and a coefficient of variation (COV) of 5. In the context of drift diffusion, L is the initial distance from the boundary, V is the drift velocity and D is the diffusion constant. (1) To evaluate the Laplace transform of the inverse Gaussian distribution numerically, simulations of first-passage times (FPT) were performed in Python using Scipy S7 . Trajectories at a reset rate r * were obtained in the following way: At each step, we sampled a new passage time τ passage and a new reset time τ reset from their corresponding distributions. If we found τ passage < τ reset , it meant that there has been a passage before the next reset time. τ passage was added to the overall simulation time and the simulation was stopped. Otherwise, it meant that the process was restarted before a transition occurred. τ reset was then added to the overall simulation time and we proceeded to sample new values of τ reset and τ passage . We continued this procedure until we encountered a successful passage, τ passage < τ reset .

Model potentials
Here we present the exact equations and parameters of the chosen model potentials. The parameters are given such that spatial distances are inÅ and potential energies are in units of 1 k B T for a temperature of 300 K.
The form of the two dimensional potential introduced by Gimondi et al. S8 is given in Eq. 4. Most of the parameters are taken as chosen there: x 1 = 2.5, x 2 = −2.5, σ 1 = 1.3, σ 2 = 1.3, y 1 = 0, y 2 = 0, λ 2 = 1. To make the right basin broader, we used a larger value of λ 1 = 2000, and a smaller coefficient for the y-coordinate harmonic spring. In order to achieve an accessible mean FPT in the absence of resetting, we slightly increased the coefficient for the x-coordinate harmonic spring, and lowered the barrier and chose A 1 = A 2 = 41.
The modified Wolfe-Quapp potential is of the form given in Eq. 5. We modified the original potential S9 by replacing the x coordinate with a rescaled x ′ = x/15, increasing the coefficients of the linear terms in both coordinates, and multiplying the resulting potential by a factor of 1.5. These modifications were done in order to achieve two remote, distinct sub-states with similar stability in the lower basin.

Derivation of Eq. 1 of the main text
Here, we will derive Eq. 1 from the main text, which connects the mean FPT at reset rates r, ⟨τ ⟩ r , to the FPT distribution at some reset rate r * , denoted by f r * (τ ). We begin with Eq. 6, derived by Reuveni, S10 which connects the FPT distribution without resetting of a random process, f (τ ), to its mean FPT with Poisson resetting at rate r, through S-5 which is identical to Eq. 1 of the main text for r * = 0. This equation holds for any distribution f (τ ), including the special case f (τ ) = f r * (τ ). Consequently, we may treat the process with reset rate r * as if it is an unbiased random process and ask what happens when adding a resetting procedure with rate r ′ to it. We rewrite the equation above for this special case in Eq. 7. The new notation emphasizes that r ′ is independent of r * and receives any values r ′ ≥ 0, as opposed to r in Eq. 1 of the main text, which only receives r ≥ r * . We signify the two independent resetting procedures by two subscripts in the left-hand side of the equation.
What is the distribution of the resulting resetting times? We combined two resetting procedures, each an individual Poisson process with reset times sampled from an exponential distribution, with a rate r * or r ′ . Due to the additive property of Poisson processes, this results in another Poisson process, with rate r = r * + r ′ . S11 Thus, the combined effect is equivalent to the introduction of a single resetting rate r = r * + r ′ , meaning ⟨τ ⟩ r * ,r ′ = ⟨τ ⟩ r * +r ′ = ⟨τ ⟩ r . Substituting this fact into Eq. 7 yields Eq. 1 of the main text,

Inference by extrapolation procedure
In this section we will describe in detail the inference procedure used in the main text to obtain the unbiased mean FPT values. We will also present alternative procedures we tested and compare their results. All procedures are based on the FPT distribution with Poisson resetting at reset rate r * , obtained from simulations.
The chosen procedure (which we denote as method A) begins with predicting the mean FPT at eight equally spaced rates in the vicinity of r * , ⟨τ ⟩ r * +i∆r , i = 1, 2, ..., 8, which was S-6 done using Eq. 8 (Eq. 1 of the main text). The results given in the main text use a spacing ∆r = 0.4r * between adjacent points. Next, the first four derivatives are evaluated at r = r * using a forward finite difference method, through where n is the order of the derivative and C n,i are coefficients given in table S1. S12 These derivatives are used to obtain an approximate fourth-order Taylor expansion of ⟨τ ⟩ r around The estimated unbiased mean FPT ⟨τ ⟩ 0 is simply the value of this function at r = 0.  We examined the sensitivity of the extrapolation to the selection of ∆r, and found less sensitivity than in other methods. Fig.S1 (a) shows the predicted ⟨τ ⟩ 0 against 1/r * for different selected values of ∆r, for the inverse Gaussian distribution with ⟨τ ⟩ 0 = 1000 ps using the analytical Laplace transform in Eq. 8. It demonstrates that ∆r values of different orders of magnitude yield similar predictions.
We also examined fitting directly the mean FPT values at different rates. We will refer to this method as method B. Here, we used a Padé approximant of the form ⟨τ ⟩ r = ar 3 +br 2 +cr+d er 2 +f r+1 , which diverges in the limit r → ∞ as does ⟨τ ⟩ r . We fitted the function numerically using Scipy S7 and substituted r = 0 to obtain the unbiased mean FPT.
A third method, which we denote as method C, uses a known connection between the mean of a distribution and its Laplace transform, ⟨τ ⟩ = − df (s) ds s=0 . We rewrite Eq. 6 as S-7 f (r) = 1 1+r⟨τ ⟩r and obtainf (r) for several rates r > r * using ⟨τ ⟩ r>r * . Then, we fit numerically a Padé approximant to these selected values of the Laplace tranform. We use the functioñ f (s) = a+bs 2 +cs a+ds 3 +es 2 +f s , which fulfills two general properties of Laplace transforms,f (0) = 1 and lim s→∞f (s) = 0. Finally, we evaluate the unbiased mean FPT using the derivative of the fitted function at zero, ⟨τ ⟩ 0 = f −c a . Fig. S1 (b) compares between the methods. Though methods B and C do not require equal spacing, here we used eight equally spaced points as needed for method A. Additional tests did not show any better performance for different spacings or number of points. We used ∆r = r * since methods B and C proved to be more sensitive to the selection of ∆r, and performed well for this value. These methods gave similar predictions to those of method A for most values of r * , but deviated for others. Since the predictions of method A improved systematically as r * → 0, as opposed to the predictions of methods B and C, we choose to present this method in the main text. Figure S1: (a) Predictions of ⟨τ ⟩ r=0 against 1/r * using method A, for different values of ∆r. (b) Predictions of ⟨τ ⟩ r=0 against 1/r * using method A, B and C. The predictions were made using exact values ⟨τ ⟩ r * +i∆r for the inverse Gaussian distribution.  The results of the main text used fixed spatial initial conditions. This is equivalent to sampling the positions initially, and after each reset, from a delta function distribution. We  The mean FPT is greater than the one achieved with fixed initial positions, and the COV is lower. As expected, the speedups are lower as well, because there is a significant probability to initiate the simulations very far from the barrier. Nevertheless, The COV remained greater than one and speedup was obtained using stochastic resetting in all model systems. This verifies that the acceleration gained by SR is not dependent on using a single specific initial condition.